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1. Introduction 

Herwig++ [1] is a general purpose Monte Carlo event generator used for simulating hard 
lepton-lepton, lepton-hadron and hadron-hadron collisions. It uses the parton shower ap- 
proach for initial and final state parton branching processes including colour coherence 
effects and azimuthal correlations. One example of a process modelled by Herwig++ is 
e^e~ annihilation to qq to form two jets (Figure 1). The jet topology (the number of jets) 
is determined by the hard cross-section of the process whilst the jet structure is deter- 
mined by Herwig++ by simulating soft and collinear branching from the primary partons 
and the conversion of the partonic final states into hadrons (hadronization) . Another 




Figure 1: 2 jet formation for e+e annihilation 
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Figure 2: LO diagram for Drell-Yan lepton pair production. 



process modelled by Herwig++ is Drell-Yan lepton pair production from hadron-hadron 
collisions which is illustrated at leading order in Figure 2. Different methods of matching 
next-to-leading order calculations to parton shower generators have been proposed and 
implemented [2, 3, 4, 5, 6, 7]. The aim of this paper is to extend the parton shower sim- 
ulation to next-to-leading order using the MC@NLO method to include the formation of 
an extra jet and NLO virtual corrections without any double counting of events. This is 
illustrated for e+e~ annihilation in Figure 3. For Drell-Yan lepton pair production, there 



are 2 real emission contributions at next-to-leading order. They are the emission of a gluon, 
q-\- q^V ^ g and the QCD Compton subprocess, q-\- g -\- q. Both are illustrated in 
Figure 4. 



e' 




Figure 3: 3 jet formation 







Figure 4: NLO diagrams for lepton pair production. 
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The generic MC@NLO method is described in [8] and has previously been successfully 
applied to the hadroproduction of gauge boson pairs [8, 9], heavy quark-antiquark pairs 
[10] and single-top production [11]. In these apphcations, the Fortran Monte Carlo event 
generator HERWIG [12] was used to simulate the parton shower and hadronization. In this 
paper however, the MC@NLO method is applied to the e+e" annihilation and Drell-Yan 
processes using Herwig++, a redeveloped version which implements new shower variables 
and an improved hadronization model [13]. 



2. e'^e annihilation 

In the massless limit, the 3-particle cross-section for the process, e'^e~ — 7* 
in Figure 5 



qqg shown 




Figure 5: Feynman diagrams for e'^e — > qqg 
is given by (neglecting Z boson exchange contributions for the moment) 

(7"-=' = ao I aXqGX^ 

where 



a^^S = ao dXqdXq — CFM{Xg,Xq) 



M{Xq,Xq) 



xl + xl 



{1 - Xq){l ~ Xq) 
2En 



2 

" 3s 



(2.1) 



(2.2) 



Cp = 4/3 (for the case of SU(3)colour representations), and the integration region is: 

< Xq, Xq < l,Xq + Xq > 1 [14]. 
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The integrand in (2.1) is divergent at Xq,Xq = 1 where the gluon is colhnear with the 
quark or antiquark or where the gluon is soft. As we shall see shortly, these singularities are 
cancelled out in the total cross-section to next-to-leading order in as- Using dimensional 
regularization, (2.1) can be evaluated to give, 



a 



ma 



(e) 



2 3 19 2 N 
^ + - + --vr^ + 0(e) 



(2.3) 



where 



H{e) = 



2(4-^), 

3(1 -e) 
(3 - 2e)r(2 - 2e) 
1 + 0(e) 



(47r: 



,2e 



(2.4) 



and n = number of dimensions. 

The total cross-section to order as is obtained by adding the contributions from the 
leading order and virtual gluon Feynman diagrams in Figure 6 to (2.3). This contribution 






Figure 6: Leading order and virtual gluon Feynman diagrams 
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1 + 



Cpas 
2tt 



Hie) - - 



+ vr2+0(e) 



(2.5) 



Taking Cp = 4/3 (for the SU(3) colour group), the total cross-section is found to be [14], 



"■total = 0-0 



l + ^ + 0{as') 

TT 



(2.6) 



For massless partons, the QCD correction at 0{as) is independent of the nature of the 
exchanged boson. Hence at this order for Z boson exchange, 



"total = ""f [1 + 



TT 



(2.7) 
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where 



2,^2 



2, 47ra K 



V2GfM^ 



(2.8) 



Aq,Vq,Ae and Ve are the axial and vector coupUng constants for the quarks and Icptons 
respectively, Gp is the Fermi constant, Mz is the pole mass of the Z boson and F is the 
total decay width of the Z boson. Note that this is the cross-section at y/s = Mz- 

As was mentioned earlier, it can be seen that the coUinear and soft singularities in 
(2.3) have been cancelled in crtotai- However, writing (Ttotai as a sum of two separate parts 
(2.3) and (2.5), makes it possible to generate a Monte-Carlo events in Xq,Xq space which 
can be fed into Herwig++ to simulate 2 and 3-jet processes. This is the subject of section 3. 

3. MC@NLO method 

Writing (Ttotai explicitly as the sum of equations (2.3) and (2.5) gives 



Ctotal — ctq 

Cpas 



(3.1) 



This can be re-written in integral form as 



O"total = (^0 dXqdx. 



2 - {M{Xq, Xq, e) - 3} + ^CpMiXq, Xq, c) 

ZTT ZTT 



(3.2) 



where 

M{Xq,Xq,€) 



//(( 



[il-Xq){l-Xq){l-Xg)r 



{1-Xq){l-Xq) 



2e 



(3.3) 



is the e+e — > qqg hard matrix element and Now, if we define a functional 

Fi as the functional which represents hadronic final states generated by the parton shower 
starting from a configuration i, a generating functional for the process e+e" hadrons 
can be written as 



F = ao J dXqdXq Fqq | 



2 - ^Cf {MiXq, Xq) - 3)} + Fqqg^CFM{Xq,Xq) 



27r 



(3.4) 



where Fqq is the functional representing the shower final states resulting from the process 

e'^e^ qq and Fg,-jg represents the final states from e"'"e~ qqg. We have set e = so 
that H(e) = 1 and M{xq,Xq,€) = M{xq,Xq) which is defined in (2.1). 

This would be wrong however because configurations starting with qq, would also 
radiate quasi-collinear gluons, with a distribution, Mcixq, Xq) given by the parton shower. 
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Likewise, configurations with qqg would generate qq-like configurations if the gluon is quasi- 
colhnear to the quark or antiquark. Mc{xq, Xg) is the parton shower branching cross- 
section, which in the massless case is given in Herwig++ by 

-I \ f Xq+Xq — ^ \ 

-C^McK,x,) = (3.5) 

for a gluon branching quasi-collinearly off the antiquark. (Interchange Xq and Xq for a gluon 
branching off the quark). This can be derived from the quasi-collinear splitting function 
defined terms of the Herwig++ evolution variables, z and q in (3.6) [15]: 



dP{q qg) = ^as[z\l - zff]^-^[l + z' - (3.6) 
ZTT q"^ 1 — z zq"^ 



where z is the momentum fraction of the quark after gluon emission relative to the parent 
quark and q is an angular variable related to the relative transverse momentum, pT of the 
quark after gluon emission via: 

I PT 1= v^(l-^)2)(zV-/i')-^Q/ , (3.7) 

/X = max (m, Qg) and Qg^ is the minimum virtuality for the quarks and gluons which is 
required to define a resolvable emission. The Dalitz plot variables Xq and Xq are related to 
the evolution variables z and q via; 

Xq = 1 — z{l — z)k , 

Xq = {2- Xq)r + (z - r) sJxl-Ap (3.8) 

where 

2 

m 

P= — , 



2 V l + p-Xq 

Xq - (2- Xq)r 

^2 



r + 



K = ^ . (3.9) 

By changing the evolution variables in (3.6) to the Dalitz plot variables in the limit where 
m = p = 0, (3.5) can be derived. The Jacobian factor for the transformation is 



z{l-z)^xl~Ap. (3.10) 

The equations given above apply to a radiating antiquark. For a radiating quark, inter- 
change Xq and Xq. Imposing the condition 

1 



Ac<-(i + yr^) (3.11) 



-7- 



defines the regions of the phase space covered by the parton showers i.e the quark and 
antiquark jets. In the massless limit this yields the function Qps in (3.12) which defines 
the phase space regions Jg, Jq and D in Figure 7. 

QpS = e[(l - Xq){Xq+Xq - 1) - X?(l - Xq)] + e[(l - Xq){Xq+X-q - 1) - X^(l - Xq)] . (3.12) 




Figure 7: Phase space for Xq,Xq showing hard (D) and soft/coUinear {Jq,Jq) gluon emission 
regions. 

The full integration region is shown shaded in the figure above. Regions Jq and Jq 
include soft and quasi-collinear gluon emission events whilst region D includes hard and 
non-collinear gluon emission events, giving rise to additional jets. 

Having defined Mc, we can now obtain the correct overall functional by subtracting 
the integral of (3.5) from the second term in (3.4) and adding it to the first term: 



I 



Fqq\2 



{2 - '^Cf (M - Mc - 3)} + Fqqg^CpiM - Mc) 



(3.13) 



This subtraction is relevant only to regions Jq and Jq in Figure 7 which include soft and 
collinear emission events. In region D, which contains hard emission events, we simply inte- 
grate the hard matrix element, M(xq,Xq) over the region. Therefore the overall generating 
functional can be written as 



dx q dx q 



+ 



D 



dx q dx q 



Fqq { 2 



Fqq [2 - ^Cf{M - Mc - 3)} + Fqqg^CpiM - Mc) 

{2-gc^(M-3)}+F, 



(3.14) 
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where J denotes the region JqU Jq. Note that the total cross-section can be retrieved by 
making the substitution Fqq, Fqqg = 1 in (3.14) as shown. 

<Ttotai = ^0 ^ dxqdxq [2 - '^Cf{M - - 3) + ^Cf{M - Mc)] 

+ j dXqdXq\2-^CF{M-?>) + ^CFM^ . (3.15) 

The 'Hit or Miss' Monte Carlo method was used to evaluate the integrals and the events 
were generated using the importance sampling method. The algorithm is described in 
Appendix A. 



4. Heavy quarks 

So far, we have discussed the limit in which the quark and anti-quark are massless. We 
shall now discuss the case where heavy flavour quarks are produced i.e. charm quarks of 
mass 1.6 GeV and bottom quarks of mass 5 GeV. There are two ways in which this can be 
treated. Both methods are described below. 

4.1 Method 1: Using the heavy quark matrix element 

The 3-particle differential cross-section for the process e^e~ -^V ^ QQd, where V repre- 
sents a vector current like the photon is given by [16]; 



1 d'^av OLs Cf 



ay dxqdxq 27r v 
where 



{xl + 2pf + (x| + 2pf + Cv 
{1 + 2p){1-xq){1-Xq) 



2p 



2p 



(1 - XQf (1 - XqY 



P 

V 



Q 



Cv = -8p(l + 2p) , 
av = cro(l -I- 2p)v . 

For an axial current contribution e^e~ A ^ QQg, we have 



1 d'^aA as Cf 
a A dxqdxQ 2n v 



{xl + 2pf + (x| + 2pf + U 
v^{1-xq){1-xq) 



2p 



2p 



{l-XQf {l-XQ? 



where 



(4.1) 



(4.2) 



(4.3) 



Ca = 2p[(3 + xgf - 19 + 4p] , 
a A = cfqv^ . 



(4.4) 



ay and a a are the Born cross-sections for heavy quark production by a vector and axial 
current respectively whilst ctq is the massless quark Born cross-section. 
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Since the partons are massive, the phase space available for gluon emission is reduced. 
It is determined by the triangle relation: 



A(4-p,x|-p,x^)<0 (4.5) 
where A(a, b, c) = a? + b^+(P — 2ab—2ac—2bc. This is equivalent to satisfying the condition 

{1 - Xq){1 - Xq){xQ + XQ - 1) > p{2 - XQ - Xgf (4.6) 

in the phase space. 

Just as in the massless limit, the phase space region can again be split into 2 regions Jq 
and Jq, containing soft and quasi-collinear gluon emission events and a region D containing 
hard and non-collinear emission events as shown in Figure 8. There is also an additional 
region labeled O outside the phase space which as we shall see in (4.9) generates 2-jet 
events. 

1 

















1 




\ 1 





Figure 8: Heavy flavour phase space for xq,Xq showing hard (D) and soft ( Jg, Jq) gluon emission 
regions. Not to scale 



As in the massless case, the total 3-particle cross-section to 0{as) is calculated by 
adding leading order and virtual gluon contributions to the integrals of (4.1) and (4.3) over 
the phase space in Figure 8. This yields for photon exchange [14] 



O'total = cry 



1 + Cl — 

TT . 



(4.7) 



where ci ~ 1 + \2p [17]. 

For Z boson exchange, to 0{as), the cross-section is given by [14] 



crtotai = cry 



l + Cl 



as 



TT J 



+ 0-A 



1 + d: 



as 



TT 



(4.8) 
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where di ^ I - 22p [17]. 

Following the same procedure as for the masslcss case, Monte Carlo events can be 
generated by writing cxtotai explicitly. For example, the jet generating functional, F can 
now be written as; 



F = av 
+ 



dxQdxQFgg ^2 + 3ci^Cf^ 
dxQdxQ [Fqq {2 - ^Cf{M - 3ci)} + FQQg^CFM 
+ ^ dxQdxQ [Fqq [2 - ^Cf{M -Mc- 3ci)} + FQQg^CF{M - Mc) 



(4.9) 



where J = JqU Jq, M is the differential cross-section defined in (4.1) and Mc is the heavy 
quark quasi-collinear branching probability given in (3.6). 

Just as in the massless limit, setting Fqq, FqQg = 1 in (4.9) recovers the vector exchange 
part of (Ttotai in (4.8). The only difference is the integral over the region O outside the 
phase space which is required to recover the full cross-section. As before, the coefficients 
of Fqq and Fqq^ generate 2-jet and 3-jet events respectively. Details of the evaluation 
of the integrals can be found in Appendix C. Event generation follows the same lines as 
described on Appendix A. 

4.2 Method 2: Using the massless quark matrix element 

Since the masses of the charm and bottom quarks (1.6 and 5 GeV respectively) are small 
compared to the center of mass energy at the Large Electron- Positron (LEP) collider (91.2 
GeV), they can in the first approximation be assumed to be massless. Hence, the massless 
matrix element can be used to obtain the 4-momentum distributions for charm and bottom 
quarks as described above for up, down and strange quarks. This is less rigorous than the 
method described in section 4.1 but it has the advantage of a smoother distribution of 
events due to the unweighting procedure being more efficient. This is the method used in 
this paper. 



5. Results on e+e annihilation 



The methods described above were used to generate e+e" events for comparison with LEP 1 
data. Details of the assignment of partonic final-state properties are described in Appendix 
D. Figures 9-13 show comparisons of event shape distributions obtained from our results 
and LEP 1 data. The massless quark matrix element method described in section 4.2 was 
used for heavy quark generation. Also compared are event shapes obtained from Herwig++ 
with the matrix element correction switched on. This is the method whereby emissions 
are only accepted into the dead region D of the phase space at a rate given by the matrix 
element. In both cases Herwig++ version 2.0.1 was used. The hadronization scale which is 
the scale at which the shower is turned off was set to the default value of 0.631 GeV and 
the 2-loop as value was used. 
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Figures 14-16 show comparisons of identified particle spectra from events of different 
flavour with SLD data [18]. In general we are able to give a good description of the data 
with the MC@NLO method. The MC@NLO results for the LEP event shape distributions 
do not differ greatly from the matrix element correction results. For example, the thrust 
distribution appears to suffer from the same problem of generating too much transverse 
structure, leading to less two-jet like event shapes. 

However, despite the similarity in results, we can be confident that the MC@NLO 
results are normalised to the full NLO cross-section including virtual corrections unlike the 
matrix element correction. 

The identified particle spectra includes hadron momenta distributions from heavy 
quark production. Although the results are similar, in some cases the MC@NLO results 
are slightly better than the matrix element correction results. This can be attributed to 
the better treatment of the heavy quark production cross-sections in the relevant plots. 



- 12 - 




Figure 11: C Parameter and D Parameter distributions and the high, Mhigh and low, Afiow 
hemisphere masses. Data from [19]. 
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The scaled momentum of charged particles for all events The scaled momentum of charged particles for all light quark events 




0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 



Figure 14: The scaled momentum distributions of all charged particles. 
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The scaled momentum of kaons for all events The scaled momentum of kaons for light quark events 




0.0 0.2 0.4 0.6 0.6 1.0 O.C 0.2 0.4 0.6 O.B 1.0 

Xp Xp 

The scaled momentum of kaons for charm events The scaled momentum of kaons for bottom quark events 




Figure 15: The scaled momentum distributions of kaons. 
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6. Drell-Yan lepton pair production 





6.1 Kinematics 

The method described above was then appHed to Drell-Yan electron pair production at the 
Tevatron. At leading order, the relevant subprocess is illustrated in Figure 2 is q + q ^ V. 
The invariant mass Q and the rapidity Y of the boson V can be written in terms of the 
momentum fractions of the incoming partons, Xq and Xq clS 

Q — XqXqS , 

Y = h^— (6.1) 

2 Xq ^ ' 

where S is the proton-antiproton center-of-mass energy. Inverting this, we have 

S_ 

%e-^ . (6.2) 

Next we consider the real emission subprocesses illustrated in Figure 4. It is convenient to 
express the kinematics of the NLO diagrams in terms of Mandelstam invariants. For gluon 
emission, we define 

S = {Pq+Pqf, 
t = {Pq-Pgf, 

U= {pq- Pgf (6.3) 

with Q"^ = s + 1 + u. For the QCD Compton process, these are given by 

S = {Pq+Pgf, 

t = {p'q-Pgf, 

U={pq-Pqf. (6.4) 

We can also express the kinematics in terms of the variables x and y which are given by 

X = — ; < X < 1 , 
s 

y = cos9;-l<y <\, (6.5) 

where 9 is the scattering angle of the emitted parton in the partonic center-of-mass system. 
Using these definitions, we can show that 

= XX1X2S , 

y^lln^^-j^-^^ + ^l (6.6) 
2 x22-{l-x){l-y) ^ ' 
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where xi,X2 are now the momentum fractions of the incoming partons in the NLO sub- 
process and are given in terms of Xg, Xg by 




{l-x){l-y) 
(l-x)(l + y) ' 



I - {l-x){l + y) 

At the parton level the Born cross-section for the production of a virtual photon is given 
by: 

where Q is the electron pair invariant mass. Extending this to the hadronic level, the Born 
cross-section becomes 

" / ^^1^^'1^^0^fli/A{Xqi)fqi/B{Xqi) q]KQ^ - Xg.Xg.S) (6.9) 

i 

where /^^/^^(xg., Q^) is the distribution function for parton i in the hadron A evaluated at 
the Born scale, Q. 

The differential cross-section for real gluon emission is given by: 

d-a ^ ^ ^C.^iis + tf + is + uf] (6.10) 



aodsdt Dg{xq)Dg{xq) 2tt sHu 

where Cp = 4/3 and Dq{xi) = x\fg{x-\) etc. Note that in this and subsequent equations, 
(To is actually J^i^y ■ '^^^ PDF ratio takes account of the change of kinematics from the 
Born momentum fractions Xg, xq to xi,X2. The corresponding differential cross-section for 
the QCD Compton subprocess is given by 

'''' -M,,^ "''!'!°'!!\> gr.g[,' + t^ + 20\l (6.11) 



aodsdt Dg{Xq)Dq{Xq) 2-K sH' 

where Tp = 1/2. The shower variables, z and k for the Drcll-Yan processes are discussed 
in detail in [15]. The invariant mass and rapidity of the boson are chosen to be preserved 
in the definition of the shower variables. Also discussed is the choice of the jet regions 
(where gluon emission is soft and/or collinear with the parent parton) for the quark q, and 
antiquark, q. In terms of the shower variables for the quark jet, the Mandelstam variables 
become 

Q2[i + (1 - z)k] 



s = 



z 

^2/ 



t = -Q^(l - z)k, 

u = ~{l-z)s (6.12) 

The jet region is then defined as the area of the s — t phase space where 

k = < kg . (6.13) 
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For the antiquark jet, we have t u and thus 



su 



(6.14) 



In order to ensure that the jet regions touch without overlapping we require kq = 1/kq. In 



the discussion that follows, we choose the symmetrical choice kg 



1. The jet and 



dead regions corresponding to this choice are labeled Jg, Jq and D respectively in Figure 
17. Now the gluon emission probability off the quark in the parton shower approximation 









D 


\J- 




\ q 





Figure 17: Jet and dead regions in s-t phase space fov q + q ^ V + g. Not to scale, 
of Herwig++ is 



dzdk 

which gives a differential cross-section 



27r ^k{l-z) 



= Mc = -P'?(^i)-Pg(^2) as^^ {s + u)[s^ + {s + u) 

aodsdt "^"^ Dq{Xq)Dq{Xq) 27r 



s^tu 



(6.15) 



(6.16) 



Interchange t u for the corresponding emission cross-section off the antiquark. Note 
that the parton shower approximation in (6.16) overestimates the matrix element expression 
(6.10) and becomes exact in the coUinear and soft limit t For the Compton subprocess 
q + g ^ y + q-, the parton shower approximation is given by 



d^ 



a 



aodsdt 



Dq{xi)Dg{x2) ^Sj^^ js + U)[U^ + {S + uf 
Dq{Xq)Dq{Xq) 2lT 



sH 



(6.17) 



Interchange t u for the subprocess g + q ^ V + q. In this case there is only one jet 
region which corresponds to the emitted quark being collinear with the gluon. The same 
jet definition in (6.13) is used and the corresponding region is shown in Figure 18. Now, 
we can further re-write the above expressions for the exact differential cross-sections and 
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Figure 18: Jet and dead regions in s-t phase space iov q + g —^ V + q. Not to scale. 

the parton shower approximation in terms of variables x and y. Using these definitions, 

Q^(l-.T)(l-y) 
2x 

„ = -Q'C . (6.18) 

In these variables the differential cross-sections become 
i^,(xi)ZJ,-(a;2) as - x)^ + (1 + x)^ 

^ Z?g(xi)lj,-(X2) «g (4 + (l-y)^ + 2x(l-y^) + .x2(l + y)2)(l-^ + x(l + y)) 

27r ^ 4(l-.T)(l-y2) 

Ug(xi)i^g(x2) as^ (3 + y2)(i_:r)2_2y(l-a;2) + 2(l + :c2) 



Dg{Xg)Dg{Xg) 2Tr 4(1-7/) 

_ £'g(3;i)L'g(x2) 1 + - 2xy(l + j/) + + y)^ 

r»,(:c,)L>,-(x,-) 27r x{l - y) ' ^ ' ^ 

Interchange y •s-^- — y in Mcgg for the antiquark jet and in Mqg and M^^^ for the process 
g + q ^ V + q. The corresponding jet regions of x — y phase space are shown in Figure 19. 
As we shall see in Section 7, expressing the cross-sections in these variables makes it easier 
to carry out the MC@NLO subtractions and divergence mappings. 

The scale at which the parton distribution functions Di{xi) , Di(x2) are evaluated was 
set to 

M = (6.20) 

which in terms of the Herwig++ variables is given by -^(1 — z)'^hQ'^. This is equal to 

I k_L I, the transverse momentum of the emitted parton in the partonic center of mass 
frame. This is the same scale used in the parton shower and is also the scale at which as 
was determined. The Drell-Yan cross-section can be computed in 2 different factorization 
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Figure 19: Jet and dead regions in x-y phase space tov q + q —^ V + g. Not to scale. 



schemes: the DIS scheme and the MS scheme. For comparison, both cross-sections were 
used for the event generation in this paper. The NLO parton distribution functions were 
obtained from the CTEQ5d (DIS) and CTEQ5m (MS) PDF sets which are frozen at a 
scale of 1 GeV. For M < 1 GeV, f{x, M) was set equal to /(x, 1). 

6.2 Next-to-leading order cross-section 

In the massless Hmit, we can integrate the differential cross-section (6.10) for the real 
emission process q + q — > V + g over y using the dimensional regularization scheme to 
regulate the divergences. The result is [20] 



(T 



NLO 

R 



r(i 



27r -^V/^V r(l-2e) 
2, /'ln(l - x) 



-^0(1 -X) - 
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— 2 mx 



1 — X 1 — X 

where e is as defined in (2.4) and the plus-prescription is defined by 



(6.21) 



/ dx{a{x))+b{x) = I a{x)[h{x)-h{l)]. 
Jo Jo 



(6.22) 



Now in addition to the real emission diagrams in Figure 4, we have virtual gluon correc- 
tions arising from self-energy and vertex corrections. These arc illustrated in Figure 20. 
Integrating these diagrams using dimensional regularization we obtain 



a 



NLO 

V 

<70 



r(i 



'r(l-2e) 



3 

e 



(6.23) 



The product of an infrared and collinear divergence is contained in the 1/e^ term and is 
cancelled out between the real and virtual diagrams to leave a pure collinear divergence, 
proportional to 1/e. This is cancelled out by QCD corrections to the quark distribution 
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Figure 20: Virtual gluon corrections to the quark-antiquark annihilation Born term q + q ^ V 



functions due to real and virtual gluon emission. In the DIS factorization scheme, the 
NLO distribution function evaluated at a scale //f is given in terms of the bare distribution 
function fq{x) as 



+ {l + z 



2^ fH^ 



1 + ^' ^^xn ^ 



lifc^ + ln^ 
e r(l - 2e) ^ 47r 



1 - z 



1 



+ 



2{l-z)+ 1-z 



Inz 



+ S + 2z-[l + ^]6{l-z) 



(6.24) 



where z = x/w. Combining these corrections gives the full NLO cross-section ratio [20] 
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(6.25) 



where signifies a sum over parton flavours q. The residual collinear divergence has 

2 

been cancelled as expected. The last term in the 0{as) term proportional to In ^ can be 
eliminated via the DGLAP equation [21, 22, 23, 24] which describes how the distribution 
functions evolve with the scale 
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d In /i^ '^'^ 
where the splitting function 
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(6.26) 



(6.27) 



describes the probability of a quark coming from the sphtting, q qg. We can use the 
above expressions to replace Dq{x,iJ?p) with Dq{x,iJ,^) in (6.25). The logarithmic term is 
then cancelled to give 
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(6.28) 



7. MC@NLO method 



Now by writing the virtual and PDF corrections in terms of the hard matrix element (6.10), 
we can rewrite (6.28) in integral from as 
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(7.1) 



The first term in the curly brackets is the sum of the Born term, virtual and QCD PDF 
corrections expressed as an integral over the variables x and y. Since the area of the x — y 
phase space is 2, there is a factor of 1/2 in the integrand. The remaining term is the 
real emission contribution to the cross-section. Now we can define a functional F as in 
(3.4) which represents final states generated from the 2 different starting configurations; 
q + q ^ V and q + q^V + gas 
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(7.2) 



where Fy and Fyg are functionals which represent final states generated from q + q ^ V 
and q + q ^ V + g starting configurations respectively. As discussed in Section 3, this is 
not entirely correct because of double counting in the final states represented by Fy arising 
from the parton shower. To resolve this issue, we subtract the contribution from the parton 
shower contributions, M^^- from the integrals in jet regions Jq and Jq in Figure 19 and 
integrate the full matrix element Mqq over the hard emission region, D. The modified 
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generating functional then becomes 
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(7.3) 



where J = JgU Jg. A similar procedure can be adopted for the Compton subprocess which 
as discussed has one jet region. In this case the QCD PDF corrections cancel out the 
collinear divergence in the matrix element. The final result is 
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(7.4) 



where Tp = 1/2 and Fyq is the functional which represents final states generated from a 
q + g ^ V + q starting configuration. Details of the algorithm used for event generation 
can be found in Appendix F. 



8. Intrinsic Pt 

In QCD, the transverse momentum of partons arises in two ways. The first which has 
been discussed above is due to the real emission of gluons and involves large momentum 
transfers. This is often termed the perturbative component of the transverse momentum 
and at large px behaves as l/px"^- At low px values, with the resummation of the double 
logarithms from soft gluon emission (as is done in parton shower generators like Herwig++) , 
this component of the pT vanishes as pT 0. 

The second way in which partons acquire transverse momentum is non-perturbative. 
It involves small momentum transfers and cannot be calculated by perturbation theory. 
Therefore this has to be modelled to fit the observed data at low pT values. A small part of 
this contribution can be attributed to quarks being confined in the transverse direction to 
within the radius of the proton and therefore gaining some intrinsic transverse momentum 
due to the uncertainty principle. 
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Data on lepton pair transverse momentum from the CFS collaboration [25] suggest 
that a Gaussian distribution (8.1) best describes the intrinsic pT distribution. 

h{pT) = -^—e^'^) (8.1) 

where PT^ms root-mean-square pT of the Gaussian distribution. 

The intrinsic pT (8.1) was implemented in Herwig++. The intrinsic pT component is 
generated according to the distribution using a random number generator and added to 
the parent partons obtained at the end of the space-like shower originating from the hard 
QCD process. 

Figure 21 shows the distribution of the per degree of freedom obtained for different 
PTrnis fits to run I (1800 GeV) CDF data [26] for Drell-Yan Z boson production and the 
best fit value can be seen to be PT^ma = 2.1 GeV. 



Chi-squareds for pTrms values for Zboson production 



5 - 




2,0 2,5 

pTrms 



3.0 



Figure 21: per degree of freedom. 



9. Results on Drell-Yan production 



Details of event generation and partonic final state properties are described in Appendices 
E, F and G. Once generated the events were showered using Herwig++ version 2.0.1 and the 
distribution of the transverse momentum of the Z boson was obtained. The hadronization 
scale was set to the default value of 0.631 GeV and the 2-loop as value was used. Figure 
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pT of the Z boson 




Figure 22: Transverse momentum of the Z boson. Data from [26]. 

22 shows a comparison of the distributions obtained for both factorization schemes with 
CDF Run I data [26]. 

Also shown in Figure 23 are the rapidity distributions of the Z boson and the positively 
charged lepton arising from its decay, for comparison distributions are compared against 
Run II {^/s = 1.96 TeV) data from the DO collaboration [27]. 

As can be seen in Figure 22 the MC@NLO method provides a good description of 
the CDF data for the transverse momentum of the Z boson. It also proves to be stable 
with respect to change of scheme. Figure 22 also shows the effect of the pT distribution 
in the low px region. The red dashed line corresponds to setting PTrms — GeV whilst 
the black line corresponds to setting pTrms = 2.1 GeV. Comparing the two, one can see the 
effect of adding the non-perturbative intrinsic pT to the parton shower which gives a better 
description of the data. 

In addition, Figure 23 shows that the predicted rapidity distributions of the Z boson 
and the positive lepton produced are stable with respect to the change of scheme. 

10. Summary and conclusions 

We have successfully applied the MC@NLO method to e^e~ annihilation and Drell-Yan 
processes modelled by the Herwig++ event generator. In general, we conclude that the 
MC@NLO method provides an improved description of event shape distributions when 
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compared to pure leading order Monte Carlo results. As we have seen, the MC@NLO 
results for e~^e~ annihilation do not differ greatly from the matrix element correction 
results although we are now confident that the results are normalised to the full NLO 
cross-section. Better results were obtained using the positive weight Nason@NLO method 
where the hardest gluon emission was generated first [5]. The MC@NLO method applied 
to the Drell-Yan process gives a good description of the transverse momentum and stable 
predictions for Z boson and lepton. The differences between the performance of the method 
for e~^e~ annihilation and Drell-Yan production may be attributed to the different shower 
variables for initial state and final state radiation. 

In addition, we have also successfully implemented the intrinsic pT component for 
hadron-hadron collisions into Herwig++. 
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A. Monte Carlo algorithm for e+e annihilation 

The integrals in (3.15) can be evaluated using a variety of Monte Carlo methods. In this 
report, the 'Hit or Miss' Monte Carlo method is used. This is the simplest and oldest form 
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of Monte Carlo integration and essentially involves finding the area of a region in phase 
space by integrating over a larger region, a binary function which is 1 in the region and 
elsewhere. The sampling method used for the points Xq,Xq is the importance sampling 
method whereby more samples are taken from regions where the integrand is large and 
less from regions where it is small. This ensures that the sampled points have the same 
distribution as the integrand. 

First, let us investigate what each of the terms in (3.15) signify. The program Herwig++ 
generates n-jets from an inputed set of n momentum space points (xj, . . . So for 2 

and 3-jet formation, we require sets of 2 and 3 momentum space points with each set 
corresponding to an event in the phase space. The relative numbers of these events as well 
as the Xi values can be obtained from (3.15) as outlined. (2-jet events have Xq,Xq = 1). In 
the discussion that follows, -y/s was set to Mz = 91.2 GeV and as{Mz) = 0.118. 

1. Randomly sample points Xq, Xq, in each of regions Jq, Jq and D of the phase space and 

(2) (3) (2) 

using the "Hit Or Miss" Monte Carlo method, evaluate the 4 integrals, I j , I j , Ij^ 
and I^' as well as their absolute sum, /. 
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r{3) 1 



(A.l) 



(3) (3) 

Note also the maximum values of the integrands in Ij and . 



2. The eventual proportion of 2-jet Monte Carlo events will be determined by the ratio 

|j(2)| + |j(2)| 

j- ° . Likewise, the proportions of 3-jet events in the soft regions Jq, Jq and the 

|j(3)| |j(3)| 

hard region D are determined by the ratios '-^j-^ and '-^-^ respectively. The algorithm 
below is then used to importance-sample the 3-jet events so that the corresponding 
{xq, Xq) values of the Monte Carlo events have the same distribution as the integrands 
in if and /g^ 

(a) For event generation in region R {R = D,Jq or Jq), randomly select a point 
Xq,Xq in that region. 



(b) Evaluate the absolute value of the integrand in for this point, | w{xq,Xq) \. 
Is I w{xq,Xq) | > i? | it;inax | ? (i? is a random number between and 1 and 
I Wraax \ IS the maximum value of | w{xq,Xq) \ determined in Step 1). 
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(c) If NO, return to (a). If YES, accept the event and set ■w^''^™= sgn w{xq,Xq) 
i.e. w^^™ = 1 if w{xq^Xq) is positive and —1 if negative. (In regions Jq and 

(3) 

Jq, M < Mc, hence the integrands and the integral, Ij in these regions are 
negative). This process is called unweighting. 

(d) Repeat the process until the correct proportion of 2-jet and 3-jet events have 
been generated. 

(2 3) 

(e) Using the importance-sampled points, obtain an estimate for the integral, /)j = 

^unw 

^^-j^ — X /, where N is the total number of Monte Carlo events generated. We 
typically use iV = 10^. 



This method of generating Monte Carlo events is termed 'Monte Carlo at Next to Leading 
Order' or MC@NLO. In this way, for a total of N events, the correct proportion of 2-jet and 
3-jet events with zb unit weight is generated with the same distribution as the integrands 
in (A.l). All of these integrals are finite, but the integrands are divergent at isolated points 
within the integration regions. Before the sampling could be carried out, the divergences 
in the integrands (which cause problems in the sampling process) had to be taken care of. 
This is the described in section B. 

B. Divergences and mappings for e+e" annihilation 



B.l Divergences in dead region D 

In region D, the hard matrix element: 

M{Xq,Xq) 



X" 



{l-Xq){l-Xq) 



(B.l) 



diverges as {xq,Xq) — > (1,0), (0,1) and (1,1). To avoid these divergences, one can map 
the divergent regions into another region in such a way that the divergence is regularized. 
This is ensured by the fact that the region of integration vanishes as the singularity is 
approached. The mappings used are presented. 



B.1.1 Region D : (1, 1) 

There is a double pole in M at {xq,Xq) 
I is mapped into a region which includes D but whose width vanishes quadratically as 



(1, 1). To avoid this pole, the region Xq,Xq > 



Xq,Xq — s- 1 [15]. The mapping used is: 



7 

4"^ 



X, 



'q = l-2{l-Xq) 



-{1-Xq 



1 3 _ 



(B.2) 
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when Xq > Xq > |. This mapping also introduces an extra weight factor of 2(1 — x'^) 
in the integrand. (Interchange Xq and Xq in both the mapping and weight factor when 
Xq > Xq > j). Figure 24 shows the region mapped (solid) and the region mapped onto 
(dashed). 




Figure 24: The mapped region (solid) and the region mapped onto (dashed) 
B.1.2 Region D : (1,0), (0,1) 

A simple pole is approached in M as {xq,Xq) — ^ (1,0) and (0, 1). In the region Xq < |, 
Xq > |, a new set of random points is generated for points which fall between the lines 
Xq = 2.5(1 — Xq) and Xq = 1 — Xq. These points have an extra weight factor to cancel the 
divergence as the pole is approached. The mapping used is: 

Xq = l- 0.25r2 , 

= (l + 1.5ri)(l-4) (B.3) 

where ri and r2 are random numbers in the range [0, 1]. This mapping introduces a weight 
factor of 2r2 in the integrand. Interchange Xq and Xq in the mapping for the region where 
Xq < |,.Tq > |. The mapped regions are shown with solid boundaries in Figure 25. 

B.2 Divergences in jet regions Jq and Jg 

In both regions Jq and Jq, there is a simple pole in the term (M — Mc) at (xq, Xq) = (1, 1). 
In the region Xq,Xq > |, a new set of random points are generated which have a weight 
factor to cancel the divergence. The mapping used in region Jq where Xq > Xq is: 

x'g = I- 0.25ri , 

4 = 1 - (1 - Xq)r2 (B.4) 
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Figure 25: Mapped regions 

where ri and r2 are random numbers in the range [0, 1]. The weight factor for this mapping 
is 2ri. For region Jg, where Xg > xq, interchange Xg and Xg in the mapping. The mapped 
regions are shown with sohd boundaries in Figure 26. There are no poles at (1,0) and (0,1) 




Figure 26: Mapped regions 



because at these points the singularities in M and Mq are cancelled out in the subtraction. 



B.3 Mapping method 

To illustrate the method behind determining the mappings, the mapping in section B.1.2 
will be explicitly calculated here. Figure 25 shows the mapped region which is between 
the lines xq = 2.5(1 — Xq) and Xq = 1 — Xq. The integrand in this region goes as and 
hence is divergent as Xq 1. By performing the change of variables shown in (B.5), this 
divergence can be regularized. 

rl r2.b{l-x^) 
I dXq / dXqf{Xq,Xq) 
Jo.75 Jl-Xq 

= [ dr2 [ dri(0.252 X 1.5r2)/(xg = l-0.25r2,Xg = (l + 1.5ri)(l-a;q-)) . (B.5) 
Jo Jo 

ri and r2 are random numbers between [0, 1] as discussed in section B.1.2 and the factor 
0.25^ X 1.5r2 is the Jacobian factor arising from the change of variables. It can be seen 
that this Jacobian factor explicitly removes the divergence in the integrand. Since 
the 'Hit Or Miss' Monte Carlo method is used to perform the integration, the weights of 
points in the mapped region must be multiplied by an area factor equal to inverse of the 
area of the mapped region i.e. ^ in this case. Hence the weight of a sampled point in this 
region is 

A 

0.25^ X 1.5r2 X — X / = 2r2/ . (B.6) 

A similar treatment is followed for all the other divergences. Now that we have taken care 
of the divergences, we need to check that the mappings give the true value of the integral. 



B.4 Testing the mappings 

The integration package VEGAS [28] was used to obtain estimates for the integrals in (3.15). 
VEGAS is an iterative and adaptive Monte Carlo integration algorithm which is based on 

importance sampling and is good for the evaluation of multidimensional integrals. 

(s) (s) 

The algorithm was used to compute integrals Ij and in section A. Forty iterations 
were carried out with 10^ program calls per iteration. The results are compared with the 
values obtained with the mappings and unweighting procedure described in sections 3 and 
B. Also worth comparing are the integrals obtained before the unweighting in step 2 of 
section A. This gives a measure of how efficient the unweighting process is. These are 
presented in the Table 1. The errors in the VEGAS results are not to be trusted as they 
appear to be underestimated. 

The estimated value for the ratio, after the mappings can also be compared with 

the true value, [1 + — ] to 0{as). This ratio is the sum of the four integrals, I^p , I^p , Ij^^ 

(3) 

and outlined in section A. Table 3 shows the relative number of 2-jet and 3-jet events 
generated from a total of 10^ events. 
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INTEGRAL 


A3) 


r(3) 


With mappings 
VEGAS 
Unweighted events 


-0.0393 ± 1 X 10-'' 
-0.03923063±6 X lO"* 
-0.0393 ± 1 X 10"^ 


0.0339 ± 1 X 10-* 
0.0339753 ± 9 x 10"^ 
0.0339 ± 1 X 10"^ 



Table 1: Comparison of VEGAS with weighted and mapping integrals 



True cross-section ratio 


Estimated cross-section ratio 


1.0375 


1.0367 ± 1.3 X 10-3 



Table 2: Comparison of the cross-sciction ratio with the true cross-section ratio to 0{as) 



Number of 2-jet events 


Number of 3- jet events 


934,567 


65,433 



Table 3: Relative number of 2-jet and 3-jet events per 1,000,000 events 
C. Heavy quark integrals 

Event generation for heavy quark production follows the same lines as discussed in Section 
A for the massless case and the soft divergence mappings used are also implemented in 
the heavy parton case. 10^ events were generated in this way to obtain a better estimate 
of the cross-section. The 3-jet integrals obtained are presented in Tables 4 and 5 for both 
charm and bottom quark production from vector boson exchange. The VEGAS results are 
also presented for comparison. The maximum absolute weights in the two regions J and 
D are also presented. 



Integral 


r(3) 


A3) 


With mappings 
VEGAS 
Maximum | weight] 


-0.0403±3 X 10--* 
-0.03989509±7x 10"^ 
2360 


0.03215±4 X 10-5 
0.0321728±7 X 10-^ 
15 


Table 4: Comparison of VEGAS and mapping integrals for charm quarks 


Integral 


r(3) 


r{3) 

-'d 


With mappings 
VEGAS 
Maximum | weight] 


-0.0483 ± 2 X 10-3 
-0.0459681 ± 6 X 10"^ 
13338 


0.02809 ±2 X 10-3 
0.0281315 ± 1 X 10-^ 
9.3 



Table 5: Comparison of VEGAS and mapping integrals for bottom quarks 

/j^'* and here are the corresponding heavy parton integrals to the integrals and 

(3) 

/n discussed in section A in the massless limit. The estimated value for the ratio, 
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after the mappings can also be compared with the true value at 0{as), 1 + ci^. 



Quark flavour 


True ratio 


Estimated ratio 


Charm 
Bottom 


1.0373 
1.0385 


1.0374±6 X 10-* 
1.0387±3 X 10-3 



Table 6: Comparison of the cross-section ratio with the true cross-section ratio to 0{as) 

The problem with this method arises during the unweighting process described in Section 
A. The efficiency of the unweighting process for a particular can be defined as the ratio of 
the integral over the region to the maximum value of its integrand — ^ . This is a measure 
of the rejection rate in the unweighting process. So the smaller the value of Wmax, the 
greater the unweighting efficiency. Now, the divergence mappings described in section B.2 
are not as effective in smoothing out the distribution of event weights in regions Jq and Jq 
(where there is a soft gluon singularity) as was the case in the massless limit. This results 
in a relatively peaked distribution in the soft gluon region and hence a relatively large 
absolute value for the maximum weight (see Tables 4 and 5). Due to this, the unweighting 
process is comparatively inefficient. 

To resolve this, we can impose a limit on the gluon 'softness' allowed in 3-jet events 
from the soft and collinear regions, Jq and Jq. 3-jet events with gluon energy fractions, 
Xg below this limit arc then considered to be 2-jet events with xq and xq equal to 1 and 
Xg = 0. For charm and bottom quarks, a limiting value of 1 x 10"^ was chosen such that 
the maximum absolute weight is sufficiently lowered (giving a smoother distribution) whilst 
the estimated cross-section is not too far off from the true cross-section (Table 7). The 
same method was applied for the axial vector boson coupling. 

D. Assigning parton properties 

Having chosen a phase space point, to generate a full event we have to assign full 4- 
momenta, as well as flavour, spin and colour information to the partons. In this section we 
address these issues in turn. 

D.l Momentum 4- vectors 

Figure 27 shows the production of a quark of 4-momentum pi , an antiquark of 4-momcntum 
P2 and a gluon of 4-momentum from an e~^e~ annihilation reaction in the centre-of-mass 
frame. It can be shown that the angles Oij between partons of momenta pi and pj satisfy 



Quark flavour 


Maximum 'soft' weight 


Estimated ratio 


Charm 
Bottom 


4.8 
43.1 


1.0381 ±5 X lO"'^ 
1.0439 ±2 X 10-3 



Table 7: Cross-section ratio and maximum weights after phase-space cut-offs 
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Figure 27: Annihilation in centre-of-mass frame 



the relations, 



cos ^12 = 
cos 013 = 
cos 6*23 = 



XqXq - 2(1 - Xg) + 4.p 



(x2-4/9)(x?-4/9) 
XqXg 2(1 Xq) 



Xg^Xq-Ap 



(D.l) 



The angular differential cross-section for the process e'^e~ — qq is defined in (D.2) [17]. 
This is the distribution of the angle 9 between the initial qq axis (before gluon emission) 
and the e'^e" axis. 



da 



dcos Q 



(1 + cos^ Q)au + 2 sin^ Gctl + 2 cos Gap ■ 



(D.2) 



In the Born approximation, 



au = l3avv + P ctaa , 
OL = - /3^)/3ayy , 

Op = 0^(TVA 



(D.3) 
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with 



-KC? 



(yVA = — [-2g/^eA/Xl(s)+4AeFe^/V/X2(s)] , 

, , ^ s(s-M|) 



(, _ M|)2 + r|M| ' 



. - Mi)2 + riMi 

(D.4) 
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Gi? is the Fermi constant, a is the electromagnetic coupling, Mz and are the mass 
and total decay width of the Z boson respectively and V/ and Af are the vector and axial 
couplings of fermion, / to the Z boson. 

By applying the unweighting procedure described in step 2 of section A to (D.2), the 
angles for each event can be distributed according to angular differential cross-section. 
Since the azimuthal angles for qq production and the q ^ qg process are isotropic, 4- 
momentum vectors for the quark, antiquark and gluon can therefore straightforwardly be 
constructed for each event. 

Initially all parton 4-momenta were generated in the massless limit. Though this is not 
essential for the parton shower, the parton 3-momenta were rescaled by a common factor 
using the Newton-Raphson iteration method to give the right parton masses (for heavy 
partons) and a gluon virtuality of 0.75 GeV. 

D.2 Flavour 

To assign flavour to the quarks, we need to investigate the flavour dependence of the total 
cross-section for e+e~ qq. Integrating (D.2) over all angles @ gives; 

af = {p + ^P{l-/3^))avv + P^aAA. (D.5) 

(Tf is the contribution of a quark of flavour / to the total cross-section. Hence the quarks 
are assigned flavours according to the relative values of af. 

D.3 Spins 

To assign spins to the partons, the matrix element for e+e~ qq is calculated for all 
possible helicity configurations and its modulus squared is used as a weight in allocating 
spins to the quarks and antiquarks. In other words, the helicity configurations are allocated 
according to their contributions to the total cross-section. It is assumed that the electron 
and positron beams are unpolarised so that the helicities of the electrons and positrons are 
assigned randomly. 
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In the masslcss limit, the chiral components are the hehcity eigenstates and so when 
a Z boson (which couples to chiral components) is exchanged, the matrix elements for 
different spin configurations are straightforward to calculate. In the heavy parton case, 
the chiral components are a mixture of helicity eigenstates so that the matrix element 
calculations become more complicated. Also, there are 4 possible helicity configurations 
in the massless limit whilst there are 8 helicity configurations in the heavy parton case 
(electrons and positrons are assumed to be massless). 

As an illustration, the matrix element for the helicity configuration, e^" g-]- is 
given in (D.6). 



nn 



(i + cose)^ 



- 2eQ5cfiXiA 



-Pc 



2 

+ c- 



+ 4 



+ Pc 



n 2 



where 



2V327ra 



(D.6) 



(D.7) 



^'r-'^'li^r right and left-handed couplings of the quarks and leptons to the 

Z boson respectively and Pcm is the centre-of-mass momentum of the quark. 

So far we have assigned the particles into helicity eigenstates. This is because the 
interface of the event record to the event generator Herwig++, requires that the electrons, 
positrons and partons are in definite helicity eigenstates. In practice however, for an 
unpolarised beam of electrons and positrons, the spins of the partons resulting from the 
annihilation reaction are superpositions of the helicity eigenstates. The spin density matrix, 
p for spin-^ particles is given by; 



which can be written as 



^(1 + ai) ^{a2-ia3) 
i(a2 + ia3) i(l-ai) 



(D.8) 



where (a) is the matrix of spin expectation values, a = (ctj,, CTy, cr^) which are the Pauli 
matrices and ^ai, and ^03 are the expectation values of the spin matrices Sz,Sx,Sy 
respectively. Since the helicity eigenstates are also eigenstates of Sz (taking the z axis 
to be along the quark momentum direction) and since |(1 + ai) and |(1 — ai) are the 
probabilities of a parton being in either state, we can assign spin states to the partons by 
distributing values of ai between and 1 according to their corresponding matrix element 
contributions. 



D.4 Gluon emission 

For the 3-jet events, the gluon can be radiated either from the quark or anti-quark. In- 
tuitively, we expect the fastest or more energetic of the two partons to be the least likely 
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to radiate the gluon i.e. most likely to keep its original direction. In the massless limit, 
it has been shown that the relative probabilities of a parton retaining its direction after 
production (not emitting the gluon) are in the ratio of their respective energies squared, 
i.e oc xf [29]. This notion is used to assign mother partons to the gluons in the 3-jet events. 
(For heavy quarks, this procedure is approximate). In addition it can be shown that the 
helicity of the gluon equals the non-emitting parton's helicity which is opposite to the spin 
of the mother parton [29]. 

D.5 Colour 




q 

Figure 28: Colour Flow 

The colour of a parton flows out of its production vertex whilst its anticolour flows 
into its production vertex. The gluon carries the anticolour of the quark colour and the 
colour corresponding to the antiquark anticolour. This is illustrated in Figure 28 for the 
specific case where a f antiquark radiates a gf gluon and becoming a g antiquark in the 
process. (The colour flow lines are essentially the same for a radiating quark). 

Colours were assigned to the partons taking into account the restrictions on the colour 
flow discussed above. Note that this colour treatment is the planar approximation, which 
is correct only to 0{j^) where Nc = 3 is the number of colours. In this approximation, 
it is always possible to draw any Feynman diagram such that no colour lines need cross as 
in the figures above. 

E. Born variables and parton flavours for Drefl-Yan production 

The Born variables and Y for each event were distributed according to the Born cross- 
section in (6.9). Hence the momentum fractions Xq and Xq given in (6.2) are generated. The 
flavour of the quarks involved in each event is also determined at the Born level according 
to the product, e^^fq{xq^, Q'^)fq{xqi-, Q"^)- Colours are assigned in the planar approximation 
which is correct to 0{j^) where Nc = 3. 
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F. Monte Carlo algorithm for Drell-Yan production 



The eight integrals hsted were evaluated by the 'Hit Or Miss' method and used to generate 
a set of unweighted events using the importance sampling method described in Section A. 



r(2) 
^Jqq 



dxdy 



Dq{Xq)Dq{Xq) 



1{6{1-x) + ^Cf 



{1-x), 



6-4x + 2{l + x'^ 



/In(l — x) 



+ ( 1 + -TT^ ) <5(1 - X, 



\ 1 — X 
^S^ = E / dxdy[Mqq-McJ, 

= ^ f dxdy [^[-^9(^l'^^)^5(^2,A^^) + q^q]l 



- Mqq + Mn 



Dq{Xq)Dq{Xq) 



1{5{1-x) + ^Cf 



{l-x), 



Q-Ax + 2{l + x')^^''^'^ 



l-x 



+ (l + -7r2)<5(l-x) 



M. 



7-(3) 
^Dqq 



/ dxdyMqq, 

j{2) = ^J dxdy r^[^9(^i'/^^)^9(^2,/x^) + g ^ g] asrrj 



Dq{Xq)Dq{Xq) 



9 



„ Tp^ s t: — 5x -I — X 
27r 2 1 2 2 



r(3) 

^Jqg 



r(2) 
^Dqg 



+ {x" + (1 + X^)) ln(l -x)}-Mqg + McJ , 

^ / dxdy[Mqg - McJ , 

j dxdy r^[^9(^i'Ai^)^9(a^2,^^) + g^5] 



Dq{Xq)Dq{Xq) 



27r ^2 12 



9 2 
ox + -X 



+ (x2 + (l + a;2))ln(l-x)}-M,g] , 



r(3) 
^Dqg 



dxdyM( 
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(F.l) 



As in the e+e~ case, the I^p integrals are negative because M < Mq- The integrals are 
all finite but there are divergences in the integrands which need to be regularized to make 
the sampling process efficient. This is the subject of Section G. So far we have discussed 
event generation in the DIS scheme. In the MS factorization scheme, the full cross-section 
ratio is 



NLO 



/ dxidx2 



x[Dq{xi,fi'^)Dq{x2, fjp') + q ^ q] 



Dq{Xq)Dq{Xq) 



5{1-x) + '^Cf 



l + x^ 
'l-x 



Inx 



+ 4(1 + 



ln(l - x) 
l-x 



+ 



+ 



+ Itt^ ) <5(1 - x) 



+ 



x[Dq{xi, ii^)Dg{x2, ix^) + q^ g\as, 



+ (x^ + (l + x^))ln 



Dq{Xq)Dq{Xq) 

(1-X)2 



27r' 



Tp <^ - + 3x - -x^ 



(F.2) 
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Hence the corresponding integrals for MS scheme event generation are 



x[Dq{xi,y?)Dq{x2, /U^) + g 1 



Dq{Xq)Dq{Xq) 



6{1-x) + ^Cf[-2 



+ 4(1 + 0;^) 



2 / ln(l - x) 



+ 



+ -8 + -7r^ U(l-x 



Mqq + Mn 



1 + x/ 
1 — a; 



Inx 



r(3) 
^Jqq 



r(2) 
^Dqq 



V] / dxdy[M, 

qJj 



x[Dq{xi,n'^)Dq{x2, //^) + g ^ g] 1 

Dq{Xq)Dq{Xq) 2 



5{1-x) + ^Cf[-2 



2tt 



l + x^ 
1-x 



Inz 



+ 4(1 + x^) (^^) ^ + (-8 + ) '^(1 - 



-M 
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^Dqq = T.l dxdyMqq, 



7-{2) 

^Jqg 



r(3) 



dxdy 



x[Dq{xi, iJ?)Dg{x2, n"^) + q g] asrp 1/1,0 _l_ 2 

Dq{Xq)Dq{Xq) 27r ^ 2 \ 2 ^ ^ 2"^ 



r(2) 



+ (x2 + (l + x2))ln 
= Y,I dxdy[Mqg - McJ , 

q J J 



dxdy 



x[Dq{xi, fj?)Dg{x2, fJ?) + q^ g] as 1/1 



Dq{Xq)Dq{Xq) 



27r 2 2 



Tft: { ^ + "ix — -X 



+ (x2 + (l + x2))ln 

^Dqg = j^dxdyMqg . 



{l-xf 
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The divergences and mappings used are discussed in Section G. 



(F.3) 



G. Divergences and mappings for Drell-Yan production 

G.l Divergences in Mqq in dead region D 

In region D, the hard matrix element:; 

, . Dq{xi)Dq{x2) as ^ y^l - x f + {1 + x f 

= Dqixq)D,ix,)2^''' (l-x)(l-y^) ^^-^^ 

diverges as (x, y) (0, 1), (0, —1) and (1, 0). The mappings used to regularize these diver- 
gences are presented. 

G.1.1 Region D : (0, 1), (0, -1) 

A simple pole is approached in Mqq as (x, y) (0, 1), (0, —1). In the region x < |, y > |, a 
new set of random points are generated for events which fall between the lines x = |(1 — y) 
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and a; = 0. The mapping used and weight factor, w for these points is: 



y =1-5^, 

X = -r2(l -y), 
w = 2ri 



(G.2) 



where ri , and r2 are random numbers in the range [0:1]. For the region x < < 
interchange y —y. The mapped regions are shown with sohd boundaries in Figure 29. 



0.6 



-0.6 



0.25 
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0.25 



-0.25 



0.25 



0.6 



Figure 29: Mapped regions 

G.1.2 Region D : (1,0) 

As a; ^ 1, there is a pole in Alqq. This divergence is regularized by applying the following 
mapping and weight factor to points which fall between the lines y = |(1 — x) and y = 
within the region x > ^,0 < y < j. 

X = 1 - -ri , 
5 



y = g?'2(i -x), 

w = 2ri . 



In the region x>|,0>?/>— |, a similar mapping is used; 



(G.3) 



X = 1 - -ri , 
5 

y = -g''2(l-a;'), 
w = 2ri . 



(G.4) 



The mapped regions are shown with solid boundaries in Figure 29. 
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G.2 Divergences in Mqg in dead region D 

In region D, the hard matrix element: 



Mqgix,y) 



Dg{xi)Dgix2)as^ (3 + y2)(l 



Dq{Xq)Dq{Xq) 2lT 



Tf- 



2?/(l -x^) + 2(1 



4(1 -y) 



(G.5) 



diverges as {x,y) (0, 1). for points witliin tlie region, x < > |, a new set of random 
points are generated for events wiiich fall between the lines a; = |(1 — y) and x = 0. The 
mapping used to regularize this divergence is: 

y =1-5-1, 

' 5 , / . 
X = -^2(1 — X ) , 



w = 2ri . 

The mapped region is shown with solid boundaries in Figure 30. 

0.25 



(G.6) 




Figure 30: Mapped regions 



G.3 Divergences in {M — Mc)qq in jet region J 

In regions Jq and Jq depicted in Figure 31, the relevant integral to be evaluated is M — Mc 
which in region Jq is given by 



(M M \ D^^i)D^as {x - l)((x - 1),/ + 2y) + (3x + 5){x + 1) 

- = Dqixq)D,ix,)2^''' MIT^^ ■ (^-'^ 



This diverges as a; ^ 0. There are no singularities at the points x = l,y = ±1 because 
the divergences at these points exactly cancel between M and Mc- Points in the region 
^ < 1)2/ > i ^'iid between the lines y = and y = 1 — 2x are mapped and re- weighted 
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according to (G.8). 



y = 1 - 2x r2, 
w = 2ri . 



Interchange y — y for the region Jg. 

0.125 

0.75 




(G. 



0.125 



Figure 31: Mapped regions 



G.4 Divergences in {M — Mc)qg in jet region J 

(M - Mc)qg is given by 

(M - Mc\g = - - 3x + 1 + 3xy{x - 1) + y^{x - l)^] 



Dq{Xq)Dq{Xq) 27r 



4x 



(G.9) 

which diverges as x — in region J. This is regularized by using the same mapping 
presented in (G.8). 
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